<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>FastGapFill tutorial</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S3 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S5 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S6 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S11 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>FastGapFill tutorial</span></h1><div  class = 'S1'><span style=' font-weight: bold;'>Author(s): Ines Thiele, Ronan M. T. Fleming, Systems Biochemistry Group, LCSB, University of Luxembourg.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Reviewer(s): </span></div><div  class = 'S1'><span>In this tutorial, we show to identify deadend metabolites and blocked reactions as well as how to fill the "gap" caused by the deadend metabolites using fastGapFill.</span></div><div  class = 'S2'><img class = "imageNode" src = "" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>The function fastGapFill [1] identifies potential missing reactions from a model using KEGG as a reference database (as provided by fastGapFill [1], default: reaction.lst). Please note that the computed solutions needs to be manually evaluated for biological relevance. Also note that not all gaps will be filled, either due to missing mapping between model metabolite abbreviations and KEGG compound ID's (defined in 'KEGG_dictionary.xls') or due to missing matches in the reference database.</span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span>Note that the fastGapFill fills all gaps at once, eliminating the issue of combinatorics, which may be associated with other gap filling tools. However, the optimal solution is defined as being the most compact (i.e., the least number of reactions to be added) solution. In contrast, other methods have been developed that make use of other criteria, such as sequence evidence. </span></div><div  class = 'S2'><img class = "imageNode" src = "" width = "612" height = "338" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>fastGapFill works on compartmentalized models. The reference database (U) will be added to each compartment, thus enabling gap filling of missing reactions in each compartment, and corresponding transport reactions connecting each metabolite in the reference database and in the reconstruction from a respective compartment with the cytosol (T).</span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span>The tutorial demonstrates the use of fastGapFill with the human metabolic reconstruction, Recon 3 [2], but it can be used with any metabolic model that is consistent with the COBRA toolbox.</span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span>Note that the solutions are not automatically categorized in this tutorial but that this needs to be done manually.</span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span>This tutorial is based on [1].</span></div><h2  class = 'S3'><span>EQUIPMENT SETUP</span></h2><div  class = 'S1'><span>If necessary, initialize the cobra toolbox:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div></div></div><div  class = 'S5'><span>For solving line</span><span>ar programming problems in FBA analysis, certain solvers are required:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >changeCobraSolver (</span><span style="color: rgb(170, 4, 249);">'gurobi'</span><span >, </span><span style="color: rgb(170, 4, 249);">'all'</span><span >, 1);</span></span></div></div></div><div  class = 'S5'><span>This tutoria</span><span>l can be ru</span><span>n </span><span>with </span><span style=' font-family: monospace;'>'glpk</span><span style=' font-family: monospace;'>'</span><span> package as linear programming solver, which does not require additional instalation and configuration. However, for the analysis of large models, such as Recon 3, it is not recommended to use </span><span style=' font-family: monospace;'>'glpk</span><span style=' font-family: monospace;'>'</span><span> but rather industrial strength solvers, such as the </span><span style=' font-family: monospace;'>'gurobi'</span><span> package. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >warning </span><span style="color: rgb(170, 4, 249);">off MATLAB:subscripting:noSubscriptsSpecified</span></span></div></div></div><h2  class = 'S3'><span>TIME REQUIREMENTS</span></h2><div  class = 'S1'><span>Hours to days, depending on model size (i.e., number of intracellular compartments) and computer performance.</span></div><h2  class = 'S6'><span>PROCEDURE</span></h2><div  class = 'S1'><span>Before proceeding with the simulations, the path for the model needs to be set up:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >pathModel = </span><span style="color: rgb(170, 4, 249);">'~/work/sbgCloud/data/models/unpublished/Recon3D_models/'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >filename = </span><span style="color: rgb(170, 4, 249);">'2017_04_28_Recon3d.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >load([pathModel, filename])</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >model = modelRecon3model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">modelRecon3model</span></span></div></div></div><div  class = 'S5'><span>In this tutorial, the used model is the generic mod</span><span>el of human metabolism, the Recon 3 [2]. </span></div><h2  class = 'S6'><span></span></h2><h2  class = 'S3'><span>Identification of deadend metabolites</span></h2><ul  class = 'S10'><li  class = 'S11'><span>Detect deadend metabolites</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >outputMets = detectDeadEnds(model)</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span>print the corresponding metabolite names</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >DeadEnds = model.mets(outputMets)</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span>These metabolites are only produced or consumed in the network and the associated reactions are blocked reactions.</span></li><li  class = 'S11'><span>Identify associated reactions:</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >[rxnList, rxnFormulaList] = findRxnsFromMets(model, DeadEnds)</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span>As you can see, these metabolites have each two reactions associated. Why are they then detected as deadend metabolites?</span></li><li  class = 'S11'><span>Let's have a look at the lower and upper bounds of these reactions:</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >model.lb(find(ismember(model.rxns,rxnList)))</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >model.ub(find(ismember(model.rxns,rxnList)))</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span>The function 'detectDeadEnds' identifies all metabolites that can be only produced or consumed in the model, aka. root deadends (or 'rootGaps'). However, in some cases it may be of interest to also identify those metabolites that are produced and consumed in the network based on topology but whose reactions cannot carry flux as they are connected to reaction(s) that involve deadend metabolites (i.e., they are connected to blocked reactions).</span></li><li  class = 'S11'><span>To identify these 'downstream' metabolites use with the option 'true':</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >[allGaps, rootGaps, downstreamGaps] = gapFind(model, </span><span style="color: rgb(170, 4, 249);">'true'</span><span >);</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span>The output include the root deadend metabolites ('rootGaps'), the downstream metabolites ('downstreamGaps') as well as a combined list of these two ('allGaps').</span></li><li  class = 'S11'><span>Note that gapFind may take some computation time depending on model size. </span></li></ul><h2  class = 'S3'><span>Identification of blocked reactions</span></h2><ul  class = 'S10'><li  class = 'S11'><span>Both the root and the downstream metabolites are part of reactions that cannot carry any flux (aka. blocked reactions), under the given network topology and simulation constraints.</span></li><li  class = 'S11'><span>Run analysis for blocked reactions. The function r</span><span>eturns </span><span>a list of blocked reactions (‘</span><span>BlockedReactions’)</span><span>.</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >BlockedReactions = findBlockedReaction(model)</span></span></div></div></div><h2  class = 'S3'><span>fastGapFill</span></h2><div  class = 'S1'><span>FastGapFill allows to set different priorities for reaction types (MetabolicRxns = internal metabolic reactions in the reference database (here, KEGG), ExchangeRxns = exchange reactions, and TransportRxns = intra- and extracellular transport reactions) using the weights. The lower the weight for a reaction type, the higher is its priority. Generally, a metabolic reaction should be prioritised in a solution over transport and exchange reactions.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >weights.MetabolicRxns = 0.1; </span><span style="color: rgb(2, 128, 9);">% Kegg metabolic reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >weights.ExchangeRxns = 0.5; </span><span style="color: rgb(2, 128, 9);">% Exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >weights.TransportRxns = 10; </span><span style="color: rgb(2, 128, 9);">% Transport reactions</span></span></div></div></div><div  class = 'S5'><span>In this example, the priority is given to the metabolic reactions follow by the exchange reactions and the transport reactions, because the internal reactions are better describe and annotated, but we not exclude the others.</span></div><div  class = 'S1'><span style=' font-weight: bold; font-style: italic;'>Hint</span><span>: The performance of algorithm is best if the weighting parameter is not 1.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Prepare the output table with statistics.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >cnt = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Model name'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Size S (original model)'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of compartments'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'List of compartments'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of blocked reactions'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of solvable blocked reactions'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Size S (flux consistent)'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Size SUX (including solvable blocked reactions)'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of added reactions (all)'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of added metabolic reactions '</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of added transport reactions '</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Number of added exchange reactions '</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Time preprocessing'</span><span >;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >Stats{cnt,1} = </span><span style="color: rgb(170, 4, 249);">'Time fastGapFill'</span><span >;cnt = cnt+1;</span></span></div></div></div><div  class = 'S5'><span>In the following, we will perform the gap filling.</span></div><div  class = 'S1'><span>Initiate parameters.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >col = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >RxnList={};</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >cnt = 1;</span></span></div></div></div><div  class = 'S5'><span>Remove constraints from exchange reactions.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >EX = strmatch(</span><span style="color: rgb(170, 4, 249);">'EX_'</span><span >,model.rxns);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >model.lb(EX)=-100;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >model.ub(EX)=100;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">EX</span></span></div></div></div><div  class = 'S5'><span> Get basic model statistics.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >Stats{cnt,i+1} = filename;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[a,b] = size(model.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >Stats{cnt,i+1} = strcat(num2str(a),</span><span style="color: rgb(170, 4, 249);">'x'</span><span >,num2str(b));cnt = cnt+1;</span></span></div></div></div><div  class = 'S5'><span> List of compartments in the model that will be considered during the gap filling.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >[tok,rem] = strtok(model.mets,</span><span style="color: rgb(170, 4, 249);">'\['</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >rem = unique(rem);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(length(rem));cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Rem = rem{1};</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 2:length(rem)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    Rem = strcat(Rem,</span><span style="color: rgb(170, 4, 249);">','</span><span >,rem{j});</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = Rem;cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">Rem tok rem</span><span >;</span></span></div></div></div><h2  class = 'S6'><span>Prepare fastGapFill.</span></h2><div  class = 'S1'><span>Here, we are going to use the function prepareFastGapFill. This function creates the new model variable (consistModel), which contains the reconstruction (S in the figure), in which the reference database (U) and the transport reactions (T) have been added to each compartment. Subsequently, all reactions that cannot carry any flux in this super-reconstruction will be removed. Note that this step is time consuming.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >tic; [consistModel,consistMatricesSUX,BlockedRxns] = prepareFastGapFill(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >tpre=toc;</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span style=' font-weight: bold;'>consistModel</span><span> </span><span>– Flux consistent subnetwork of the input model</span></li><li  class = 'S11'><span style=' font-weight: bold;'>consistMatricesSUX</span><span> </span><span>– Flux consistent</span><span> </span><span style=' font-style: italic;'>SUX</span><span> </span><span>matrix, which contains the flux consistent</span><span> </span><span style=' font-style: italic;'>S</span><span> </span><span>matrix (model), the universal database placed in all cellular compartments along with transport reactions for each metabolite from cytosol to compartment and exchange reactions for all extracellular metabolites.</span></li><li  class = 'S11'><span style=' font-weight: bold;'>BlockedRxns</span><span> </span><span>– Blocked reactions in model.</span></li></ul><div  class = 'S1'><span>Add on more statistics to the table. Here, we will report how many gaps cannot be filled at all in the starting reconstruction.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(length(BlockedRxns.allRxns));cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(length(BlockedRxns.solvableRxns));cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[a,b] = size(consistModel.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = strcat(num2str(a),</span><span style="color: rgb(170, 4, 249);">'x'</span><span >,num2str(b));cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[a,b] = size(consistMatricesSUX.S);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >Stats{cnt,i+1} = strcat(num2str(a),</span><span style="color: rgb(170, 4, 249);">'x'</span><span >,num2str(b));cnt = cnt+1;</span></span></div></div></div><h2  class = 'S6'><span>Perform fastGapsFill. </span></h2><div  class = 'S1'><span>Now our model is ready to be run by the fastGapFill function. The main aim of this function is to find the most compact set of reactions, i.e., exchange, transport, and metabolic reactions from the UX matrix to be added to the input model to close the models in the gaps.</span></div><div  class = 'S1'><span>The variable epsilon defines effectively the minimum flux through a reactions that should be achieved when adding the reactions. Please refer to [3] for more details on epsilon and its use.</span></div><div  class = 'S1'><span>Note this step is time-consuming.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >epsilon = 1e-4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >tic; [AddedRxns] = fastGapFill(consistMatricesSUX,epsilon, weights);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >tgap=toc;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(length(AddedRxns.rxns));cnt = cnt+1;</span></span></div></div></div><h2  class = 'S6'><span>Postprocessing of results</span></h2><div  class = 'S1'><span>Here, we assign  generalized functions to the output from</span><span> </span><span style=' font-style: italic;'>fastGapFill</span><span> </span><span>(e.g., ‘Metabolic reaction’, ‘Transport reaction’, ‘Exchange reaction’). It also provides basic statistics for the solution. One can use the option ‘IdentifyPW’ to compute for each solved blocked reaction a flux vector, such that the flux through this reaction is maximized while the sum of all fluxes is minimized. This analysis can be helpful to put the solution reactions into the network context.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >IdentifyPW = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[AddedRxnsExtended] = postProcessGapFillSolutions(AddedRxns,model,BlockedRxns);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">AddedRxns</span><span >;</span></span></div></div></div><ul  class = 'S10'><li  class = 'S11'><span style=' font-weight: bold;'>AddedRxnsExtended</span><span>  is a s</span><span>tructure containing the information present in </span><span style=' font-style: italic;'>AddedRxns</span><span> </span><span>along with the statistics and if desired pathways containing the flux vectors.</span></li></ul><div  class = 'S1'><span>Add further statistics about computed solutions to the Stats table.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(AddedRxnsExtended.Stats.metabolicSol);cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(AddedRxnsExtended.Stats.transportSol);cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(AddedRxnsExtended.Stats.exchangeSol);cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(tpre);cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Stats{cnt,i+1} = num2str(tgap);cnt = cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">a b</span></span></div></div></div><div  class = 'S5'><span>    Generate a reaction list</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >RxnList{1,col}=filename;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.rxns; col = col + 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >RxnList{1,col}=filename;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.rxnFormula; col = col + 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >RxnList{1,col}=filename;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.subSystem; col = col + 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">AddedRxnsExtended tgap tpre j BlockedRxns i cnt consistM*</span></span></div></div></div><h2  class = 'S6'><span>Analysis of results</span></h2><div  class = 'S1'><span>Now, we are going analyze  the candidate gap filling reactions, contained in 'RxnList'.The first column correspond to the name of the reaction in the model, the second column correspond to the reactions and the last column give the compartment of the reactions.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >RxnList</span></span></div></div></div><div  class = 'S5'><span>Now, we are going to identify the deadend metabolites, for this we are going to use the function detectDeadEnds, which returns the indices of all deadend metabolites.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >metid = detectDeadEnds(model,false);</span></span></div></div></div><div  class = 'S5'><span>Print the metabolite abbreviations and metabolite names on the screen.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >horzcat(model.mets(metid),model.metNames(metid))</span></span></div></div></div><h2  class = 'S6'><span>Trouble shooting</span></h2><div  class = 'S1'><span style=' font-weight: bold; font-style: italic;'>Error</span><span>: "Combining Remaining Genetic Information: Reference to non-existent field 'rules'."</span></div><div  class = 'S1'><span style=' font-weight: bold; font-style: italic;'>Solution</span><span>: "Your model structure does not contain the rules field. The rules field specifies which genes encode for the reactions. You could create an array of the same length as the number of reactions in your model and empty fields."</span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><span style=' font-weight: bold; font-style: italic;'>Error</span><span>: " </span><span>I have downloaded the files and I am using the gurobi5 solver. When I am trying to run one of theirs runGapFill_example to check correct setup, this is what I am experiencing..</span></div><div  class = 'S1'><span>&gt;&gt; runGapFill_example</span></div><div  class = 'S1'><span>|J|=386  Undefined function 'cplexlp' for input arguments of type 'double'."</span></div><div  class = 'S1'><span style=' font-weight: bold; font-style: italic;'>Solution</span><span>: "Fastgapfill requires ibm cplx to be installed. Academic licenses are free."</span></div><h2  class = 'S3'><span>References</span></h2><div  class = 'S1'><span>[1] fastGapFill.</span></div><div  class = 'S1'><span>[2] Brunk et al, Recon 3D, submitted.</span></div><div  class = 'S1'><span>[3] fastCore</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% FastGapFill tutorial
% *Author(s): Ines Thiele, Ronan M. T. Fleming, Systems Biochemistry Group, 
% LCSB, University of Luxembourg.*
% 
% *Reviewer(s):* 
% 
% In this tutorial, we show to identify deadend metabolites and blocked reactions 
% as well as how to fill the "gap" caused by the deadend metabolites using fastGapFill.
% 
% 
% 
% The function fastGapFill [1] identifies potential missing reactions from a 
% model using KEGG as a reference database (as provided by fastGapFill [1], default: 
% reaction.lst). Please note that the computed solutions needs to be manually 
% evaluated for biological relevance. Also note that not all gaps will be filled, 
% either due to missing mapping between model metabolite abbreviations and KEGG 
% compound ID's (defined in 'KEGG_dictionary.xls') or due to missing matches in 
% the reference database.
% 
% 
% 
% Note that the fastGapFill fills all gaps at once, eliminating the issue of 
% combinatorics, which may be associated with other gap filling tools. However, 
% the optimal solution is defined as being the most compact (i.e., the least number 
% of reactions to be added) solution. In contrast, other methods have been developed 
% that make use of other criteria, such as sequence evidence. 
% 
% 
% 
% fastGapFill works on compartmentalized models. The reference database (U) 
% will be added to each compartment, thus enabling gap filling of missing reactions 
% in each compartment, and corresponding transport reactions connecting each metabolite 
% in the reference database and in the reconstruction from a respective compartment 
% with the cytosol (T).
% 
% 
% 
% The tutorial demonstrates the use of fastGapFill with the human metabolic 
% reconstruction, Recon 3 [2], but it can be used with any metabolic model that 
% is consistent with the COBRA toolbox.
% 
% 
% 
% Note that the solutions are not automatically categorized in this tutorial 
% but that this needs to be done manually.
% 
% 
% 
% This tutorial is based on [1].
%% EQUIPMENT SETUP
% If necessary, initialize the cobra toolbox:

initCobraToolbox(false) % false, as we don't want to update
%% 
% For solving linear programming problems in FBA analysis, certain solvers are 
% required:

changeCobraSolver ('gurobi', 'all', 1);
%% 
% This tutorial can be run with |'glpk'| package as linear programming solver, 
% which does not require additional instalation and configuration. However, for 
% the analysis of large models, such as Recon 3, it is not recommended to use 
% |'glpk'| but rather industrial strength solvers, such as the |'gurobi'| package. 

warning off MATLAB:subscripting:noSubscriptsSpecified
%% TIME REQUIREMENTS
% Hours to days, depending on model size (i.e., number of intracellular compartments) 
% and computer performance.
%% PROCEDURE
% Before proceeding with the simulations, the path for the model needs to be 
% set up:

pathModel = '~/work/sbgCloud/data/models/unpublished/Recon3D_models/';
filename = '2017_04_28_Recon3d.mat';
load([pathModel, filename])
model = modelRecon3model;
clear modelRecon3model
%% 
% In this tutorial, the used model is the generic model of human metabolism, 
% the Recon 3 [2]. 
%% 
%% Identification of deadend metabolites
%% 
% * Detect deadend metabolites

outputMets = detectDeadEnds(model)
%% 
% * print the corresponding metabolite names

DeadEnds = model.mets(outputMets)
%% 
% * These metabolites are only produced or consumed in the network and the associated 
% reactions are blocked reactions.
% * Identify associated reactions:

[rxnList, rxnFormulaList] = findRxnsFromMets(model, DeadEnds)
%% 
% * As you can see, these metabolites have each two reactions associated. Why 
% are they then detected as deadend metabolites?
% * Let's have a look at the lower and upper bounds of these reactions:

model.lb(find(ismember(model.rxns,rxnList)))
model.ub(find(ismember(model.rxns,rxnList)))
%% 
% * The function 'detectDeadEnds' identifies all metabolites that can be only 
% produced or consumed in the model, aka. root deadends (or 'rootGaps'). However, 
% in some cases it may be of interest to also identify those metabolites that 
% are produced and consumed in the network based on topology but whose reactions 
% cannot carry flux as they are connected to reaction(s) that involve deadend 
% metabolites (i.e., they are connected to blocked reactions).
% * To identify these 'downstream' metabolites use with the option 'true':

[allGaps, rootGaps, downstreamGaps] = gapFind(model, 'true');
%% 
% * The output include the root deadend metabolites ('rootGaps'), the downstream 
% metabolites ('downstreamGaps') as well as a combined list of these two ('allGaps').
% * Note that gapFind may take some computation time depending on model size. 
%% Identification of blocked reactions
%% 
% * Both the root and the downstream metabolites are part of reactions that 
% cannot carry any flux (aka. blocked reactions), under the given network topology 
% and simulation constraints.
% * Run analysis for blocked reactions. The function returns a list of blocked 
% reactions (‘BlockedReactions’).

BlockedReactions = findBlockedReaction(model)
%% fastGapFill
% FastGapFill allows to set different priorities for reaction types (MetabolicRxns 
% = internal metabolic reactions in the reference database (here, KEGG), ExchangeRxns 
% = exchange reactions, and TransportRxns = intra- and extracellular transport 
% reactions) using the weights. The lower the weight for a reaction type, the 
% higher is its priority. Generally, a metabolic reaction should be prioritised 
% in a solution over transport and exchange reactions.

weights.MetabolicRxns = 0.1; % Kegg metabolic reactions
weights.ExchangeRxns = 0.5; % Exchange reactions
weights.TransportRxns = 10; % Transport reactions
%% 
% In this example, the priority is given to the metabolic reactions follow by 
% the exchange reactions and the transport reactions, because the internal reactions 
% are better describe and annotated, but we not exclude the others.
% 
% _*Hint*_: The performance of algorithm is best if the weighting parameter 
% is not 1.
% 
% *Prepare the output table with statistics.*

cnt = 1;
Stats{cnt,1} = 'Model name';cnt = cnt+1;
Stats{cnt,1} = 'Size S (original model)';cnt = cnt+1;
Stats{cnt,1} = 'Number of compartments';cnt = cnt+1;
Stats{cnt,1} = 'List of compartments';cnt = cnt+1;
Stats{cnt,1} = 'Number of blocked reactions';cnt = cnt+1;
Stats{cnt,1} = 'Number of solvable blocked reactions';cnt = cnt+1;
Stats{cnt,1} = 'Size S (flux consistent)';cnt = cnt+1;
Stats{cnt,1} = 'Size SUX (including solvable blocked reactions)';cnt = cnt+1;
Stats{cnt,1} = 'Number of added reactions (all)';cnt = cnt+1;
Stats{cnt,1} = 'Number of added metabolic reactions ';cnt = cnt+1;
Stats{cnt,1} = 'Number of added transport reactions ';cnt = cnt+1;
Stats{cnt,1} = 'Number of added exchange reactions ';cnt = cnt+1;
Stats{cnt,1} = 'Time preprocessing';cnt = cnt+1;
Stats{cnt,1} = 'Time fastGapFill';cnt = cnt+1;
%% 
% In the following, we will perform the gap filling.
% 
% Initiate parameters.

col = 1;
RxnList={};
cnt = 1;
%% 
% Remove constraints from exchange reactions.

EX = strmatch('EX_',model.rxns);
model.lb(EX)=-100;
model.ub(EX)=100;
clear EX
%% 
% Get basic model statistics.

Stats{cnt,i+1} = filename;cnt = cnt+1;
[a,b] = size(model.S);
Stats{cnt,i+1} = strcat(num2str(a),'x',num2str(b));cnt = cnt+1;
%% 
% List of compartments in the model that will be considered during the gap filling.

[tok,rem] = strtok(model.mets,'\[');
rem = unique(rem);
Stats{cnt,i+1} = num2str(length(rem));cnt = cnt+1;
Rem = rem{1};
for j = 2:length(rem)
    Rem = strcat(Rem,',',rem{j});
end
Stats{cnt,i+1} = Rem;cnt = cnt+1;
clear Rem tok rem;
%% Prepare fastGapFill.
% Here, we are going to use the function prepareFastGapFill. This function creates 
% the new model variable (consistModel), which contains the reconstruction (S 
% in the figure), in which the reference database (U) and the transport reactions 
% (T) have been added to each compartment. Subsequently, all reactions that cannot 
% carry any flux in this super-reconstruction will be removed. Note that this 
% step is time consuming.

tic; [consistModel,consistMatricesSUX,BlockedRxns] = prepareFastGapFill(model);
tpre=toc;
%% 
% * *consistModel* – Flux consistent subnetwork of the input model
% * *consistMatricesSUX* – Flux consistent _SUX_ matrix, which contains the 
% flux consistent _S_ matrix (model), the universal database placed in all cellular 
% compartments along with transport reactions for each metabolite from cytosol 
% to compartment and exchange reactions for all extracellular metabolites.
% * *BlockedRxns* – Blocked reactions in model.
%% 
% Add on more statistics to the table. Here, we will report how many gaps cannot 
% be filled at all in the starting reconstruction.

Stats{cnt,i+1} = num2str(length(BlockedRxns.allRxns));cnt = cnt+1;
Stats{cnt,i+1} = num2str(length(BlockedRxns.solvableRxns));cnt = cnt+1;
[a,b] = size(consistModel.S);
Stats{cnt,i+1} = strcat(num2str(a),'x',num2str(b));cnt = cnt+1;
[a,b] = size(consistMatricesSUX.S);
Stats{cnt,i+1} = strcat(num2str(a),'x',num2str(b));cnt = cnt+1;
%% Perform fastGapsFill. 
% Now our model is ready to be run by the fastGapFill function. The main aim 
% of this function is to find the most compact set of reactions, i.e., exchange, 
% transport, and metabolic reactions from the UX matrix to be added to the input 
% model to close the models in the gaps.
% 
% The variable epsilon defines effectively the minimum flux through a reactions 
% that should be achieved when adding the reactions. Please refer to [3] for more 
% details on epsilon and its use.
% 
% Note this step is time-consuming.

epsilon = 1e-4;
tic; [AddedRxns] = fastGapFill(consistMatricesSUX,epsilon, weights);
tgap=toc;
Stats{cnt,i+1} = num2str(length(AddedRxns.rxns));cnt = cnt+1;
%% Postprocessing of results
% Here, we assign  generalized functions to the output from _fastGapFill_ (e.g., 
% ‘Metabolic reaction’, ‘Transport reaction’, ‘Exchange reaction’). It also provides 
% basic statistics for the solution. One can use the option ‘IdentifyPW’ to compute 
% for each solved blocked reaction a flux vector, such that the flux through this 
% reaction is maximized while the sum of all fluxes is minimized. This analysis 
% can be helpful to put the solution reactions into the network context.

IdentifyPW = 0;
[AddedRxnsExtended] = postProcessGapFillSolutions(AddedRxns,model,BlockedRxns);
clear AddedRxns;
%% 
% * *AddedRxnsExtended*  is a structure containing the information present in 
% _AddedRxns_ along with the statistics and if desired pathways containing the 
% flux vectors.
%% 
% Add further statistics about computed solutions to the Stats table.

Stats{cnt,i+1} = num2str(AddedRxnsExtended.Stats.metabolicSol);cnt = cnt+1;
Stats{cnt,i+1} = num2str(AddedRxnsExtended.Stats.transportSol);cnt = cnt+1;
Stats{cnt,i+1} = num2str(AddedRxnsExtended.Stats.exchangeSol);cnt = cnt+1;

Stats{cnt,i+1} = num2str(tpre);cnt = cnt+1;
Stats{cnt,i+1} = num2str(tgap);cnt = cnt+1;
clear a b
%% 
% Generate a reaction list

RxnList{1,col}=filename;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.rxns; col = col + 1;
RxnList{1,col}=filename;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.rxnFormula; col = col + 1;
RxnList{1,col}=filename;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.subSystem; col = col + 1;

clear AddedRxnsExtended tgap tpre j BlockedRxns i cnt consistM*
%% Analysis of results
% Now, we are going analyze  the candidate gap filling reactions, contained 
% in 'RxnList'.The first column correspond to the name of the reaction in the 
% model, the second column correspond to the reactions and the last column give 
% the compartment of the reactions.

RxnList
%% 
% Now, we are going to identify the deadend metabolites, for this we are going 
% to use the function detectDeadEnds, which returns the indices of all deadend 
% metabolites.

metid = detectDeadEnds(model,false);
%% 
% Print the metabolite abbreviations and metabolite names on the screen.

horzcat(model.mets(metid),model.metNames(metid))
%% Trouble shooting
% _*Error*_: "Combining Remaining Genetic Information: Reference to non-existent 
% field 'rules'."
% 
% _*Solution*_: "Your model structure does not contain the rules field. The 
% rules field specifies which genes encode for the reactions. You could create 
% an array of the same length as the number of reactions in your model and empty 
% fields."
% 
% 
% 
% _*Error*_: " I have downloaded the files and I am using the gurobi5 solver. 
% When I am trying to run one of theirs runGapFill_example to check correct setup, 
% this is what I am experiencing..
% 
% >> runGapFill_example
% 
% |J|=386  Undefined function 'cplexlp' for input arguments of type 'double'."
% 
% _*Solution*_: "Fastgapfill requires ibm cplx to be installed. Academic licenses 
% are free."
%% References
% [1] fastGapFill.
% 
% [2] Brunk et al, Recon 3D, submitted.
% 
% [3] fastCore
##### SOURCE END #####
-->
</div></body></html>